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Abstract 

Background: Recently, network meta-analysis of survival data with a multidimensional treatment effect was 
introduced. With these models the hazard ratio is not assumed to be constant over time, thereby reducing the 
possibility of violating transitivity in indirect comparisons. However, bias is still present if there are systematic 
differences in treatment effect modifiers across comparisons. 

Methods: In this paper we present multidimensional network meta-analysis models for time-to-event data that are 
extended with covariates to explain heterogeneity and adjust for confounding bias in the synthesis of evidence 
networks of randomized controlled trials. The impact of a covariate on the treatment effect can be assumed to be 
treatment specific or constant for all treatments compared. 

Results: An illustrative example analysis is presented for a network of randomized controlled trials evaluating 
different interventions for advanced melanoma. Incorporating a covariate related to the study date resulted in 
different estimates for the hazard ratios over time than an analysis without this covariate, indicating the importance 
of adjusting for changes in contextual factors over time. 

Conclusion: Adding treatment-by-covariate interactions to multidimensional network meta-analysis models for 
published survival curves can be worthwhile to explain systematic differences across comparisons, thereby reducing 
inconsistencies and bias. An additional advantage is that heterogeneity in treatment effects can be explored. 



Background 

Healthcare decision-making requires comparisons of 
relevant competing treatments for a particular disease 
state. In the absence of trials involving a direct compari- 
son of interventions, an indirect comparison can provide 
useful evidence for the treatment effects between com- 
peting interventions [1-6]. Even when direct evidence is 
available, combining consistent direct and indirect esti- 
mates will result in more precise treatment effect esti- 
mates [2-6]. In general, if the evidence base consists of 
multiple randomized controlled trials (RCTs) with each 
trial having at least one intervention in common with 
another, this network of RCTs can be synthesized with 
meta-analysis techniques: a network meta-analysis (or 
mixed treatment comparison meta-analysis) [1-6]. 
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Traditional (network) meta-analyses for survival out- 
comes are based on hazard ratios (HR) and rely on the 
proportional hazards assumption, which is implausible if 
the hazard functions of competing interventions cross. 
Ouwens et al. (2011) and Jansen (2011) presented meth- 
ods for (network) meta-analysis of survival data using a 
multidimensional treatment effect as an alternative to 
the synthesis of the constant HRs [7,8]. The hazard 
functions of the interventions in a trial are modeled 
using known parametric survival functions (e.g. Weibull 
or Gompertz) or fractional polynomials and the differ- 
ence in the parameters are considered the multidimen- 
sional treatment effect, which are synthesized (and 
indirectly compared) across studies. In essence, with this 
approach the treatment effects are represented by mul- 
tiple parameters rather than a single parameter. 

For network meta-analysis it is important that tran- 
sitivity (i.e. the consistency assumption) holds for the 
treatment effect measure of interest [2,3,9]. Violations of 
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the proportional hazards assumption will compromise 
transitivity and can result in biased indirect and mixed 
treatment comparisons of survival outcomes [8]. By 
incorporating additional parameters for the treatment 
effect as proposed with the methods by Ouwens et al 
and Jansen, the proportional hazards assumption is 
relaxed and indirect comparisons are less likely to be 
biased [7,8]. However, violation of the proportional 
hazards assumption is not the only possible source of 
bias. Since randomization of patients does not hold 
across trials, there might be an imbalance in patient or 
study level covariates across comparisons [4-6,9]. If these 
covariates are effect-modifiers of the (time varying) HRs, 
transitivity will be violated and the network meta- 
analysis will result in biased estimates, despite modeling 
a multidimensional treatment effect. 

In this paper, the models proposed by Ouwens et al. 
and Jansen are extended with treatment-by-covariate 
interactions to adjust for confounding bias due to sys- 
tematic differences in treatment effect modifiers across 
comparisons. These models also have the advantage of 
explaining heterogeneity in the treatment effects. The 
method is illustrated with an example. 

Methods 

Multidimensional network meta-analysis models for 
survival data 

If AB trials and AC trials are comparable in terms of ef- 
fect modifiers (i.e. covariates that affect the treatment ef- 
fect) then an indirect estimate for the treatment effect of 
C versus B {d BC ) can be obtained from the estimates of 
the effect of B versus A (d AB ) and the effect of C versus 
A (d AC ): d BC = d AC - d AB ; as such, transitivity holds 
[2,3,6]. For an arm-based model with the outcome of 
treatment ^asa function of time, f x (t) where x = A, B, or 
C and t represents time, this consistency assumption 
translates into: 

(f C (t) -fB{t)) = (fdt) -/AW) - Ut) ~fA{t)) 

The hazard function is of central interest to 
summarize survival or time-to-event data. The hazard 
function describes the instantaneous event (e.g. death) 
rate at any point in time. A random effects network 
meta-analysis model for hazard rates by treatment arm 
without assuming a particular distribution for the hazard 
function over time can defined according to: 



fa) 



I {i jht & = A,B,C, if k = b 

I ft fi t + Sjbk if k 'alphabetically after' b 

8jbk~normal(dbk, 0 2 ) = normal{d A k — d Ab , a 2 ) 



hazard rates at time t with comparator treatment b. Sj b k 
reflects the study specific constant log HRs of treatment 
k relative to comparator treatment b and are drawn from 
a normal distribution with the pooled estimates 
expressed in terms of the overall reference treatment 
A: d bk =d Ak - d Ab with d AA = 0. Variance o 2 reflects the 
heterogeneity across studies. As an alternative to a non- 
parametric baseline hazard function, the development of 
the hazard rate over time can be described by parametric 
functions, such as Weibull or Gompertz. 

Two-dimensional treatment effects without covariates (Model 1) 

Ouwens et al. and Jansen introduced network meta-analysis 
models for parametric hazard functions [7,8]. A network 
meta-analysis model for the hazard rate with a two dimen- 
sional treatment effect can be defined according to: 



with t° = ln(t) 



P\jb ) \d\Ak — d\Ab 

8ojbk~norrnal(d 0 Ak ~ d 0A b, & 2 ) 



fc = A,B,C, ifk = b 
if k 'alphabetically after' b 



where hj kt reflect the underlying hazard rate in trial / for 
intervention k at time point t. fy bt are the study ; specific 



where hjkt again reflects the underlying hazard rate in 
trial j for intervention k at time point t and is now 
described as a function of time t with p={-2,-l, 
-0.5,0,0.5, l,2,3j and £°=ln(t) with treatment and study 
specific scale and shape parameters fi 0 j k and fiy k . (In the 
example presented below details on the likelihood and 
data structure are provided). If fiy k equals 0, a constant 
log hazard function is obtained, reflecting exponentially 
distributed survival times. If fiy k ^ 0 and /? = 1 a linear 
log hazard function is obtained which corresponds to a 
Gompertz survival function [8]. If fiy k 7^ 0 and p = 0 
a Weibull hazard function is obtained. The vectors 

^°i b ) are trial specific and reflect the true underlying 

Pyb J 

scale and shape parameters of the comparator treat- 
ment b. S 0 j bk is the study specific difference in the 
scale parameter /? 0 of the log hazard curve for treat- 
ment k relative to comparator treatment b. S 0 j bk are 
drawn from a normal distribution with the pooled 
estimates expressed in terms of the overall reference 
treatment A: d 0Ak -d 0Ab with d 0AA = 0. The parameters 
doAk correspond to the treatment effect of k relative to 
overall reference treatment A with a proportional hazard 
model. The pooled difference in the shape parameter fi x 
of the log hazard curve for treatment k relative to com- 
parator treatment b is expressed as d 1Ak - d 1Ab with d XAA 
= 0. d\ Ak reflects the change in the log HR over 
time. For a proportional hazards model d 1Ak equals 0. By 
incorporating d 1Ak in addition to d 0Ak a multidimen- 
sional treatment effect is used. For additional flexibility, 
the first order fractional polynomial model can be 



Jansen and Cope BMC Medical Research Methodology 2012, 12:152 
http://www.biomedcentral.eom/1 471-2288/1 2/1 52 



Page 3 of 16 



generalized to a 2 n order fractional polynomial model, 
representing 3-dimensional treatment effects [8]. 

Variance Oq reflects the heterogeneity in the difference 
in the scale parameters across studies. A random effects 
model with only a heterogeneity parameter for d 0Ak 
implies that the between study variance of the log hazard 
ratios remains constant over time. 

Two-dimensional treatment effects with treatment specific 
covariate interactions (Model 2) 

Since randomization only holds within a trial, the 
distribution of covariates may vary across studies for 
a particular type of comparison, as well as between 
different types of comparisons. Variation in treatment 
effect modifiers across studies within comparisons results 
in heterogeneity and an imbalance in effect modifiers be- 
tween different types of comparisons results in violations of 
the consistency assumption. [2-4,6,9,10]. Network meta- 
analysis models can be extended with treatment-by- 
covariate interactions in an attempt to adjust for this con- 
founding bias or to explore sources of heterogeneity [9-14]. 

In line with previous meta-regression models for 
network meta-analysis where the treatment effect 
acts on one parameter [10], model 1 can be extended 
with study level covariates to explore treatment-by- 
covariate interactions. The specification assumes that 
all treatment-by-covariate interactions are different 
for each treatment k relative to overall reference 
treatment A: 

( Pojk\ = )\Pyb J & = A,B,C, \ik = b 

\Pijk J \( Ptyb\ f $ojbk \ if k ' alphabetically after' b 
J \d\Ak — diAb J 
S 0 jbk-normal(d obk + P xbk Xj,a 0 2 ) 

= normal(d 0Ak - d 0A b + {P xAk ~ PxAb) x j> °o 2 ) 

Pxbk reflects the impact of study level covariate Xj on 
the difference in the scale parameters of the hazard 
functions with treatment k relative to control treatment 
b. Now dobk is the difference in scale treatment k relative 
to b when the covariate value equals zero. Since fi xbk = 
PxAk ~ PxAb with ft xAA = 0, k - 1 different and independ- 
ent regression coefficients for f$ xAk will be estimated with 
the model. As an alternative to independent treatment- 
by-covariate interactions, one can also assume exchange- 
able interaction effects [10]. 

Two-dimensional treatment effects with constant covariate 
interaction (Model 3) 

One can also assume that the impact of the covariate 
Xj on the scale parameter of each treatment k relative 
to A is the same for all treatments. This assumption 
can be defended when treatments indirectly compared 



are all from the same class and there is no (biological) 
reason to assume that a patient characteristic, or any 
other contextual aspect of the study, modifies treat- 
ment effects differently for the different drugs com- 
pared. Furthermore, the assumption of constant 
treatment-by-covariate interaction can also be useful 
when evaluating the impact of study (design) charac- 
teristics (or bias) on treatment effects. [10,14]. The 
corresponding network meta-analysis model will be: 

ln(V) = Pojk + Pyk* with *° = ln (t) 

( Pojk\ = )\f*ijb J fc = A,B,C, if k = b 

\Pijk J \f ftojb \ _j_ f $ojbk \ if k 'alphabetically after' b 
\\f*ljb J \d\Ak — d\Ab J 

^ f normal (d 0 Ak - d 0A b + fi x Xj, a 0 2 ) if b = A 
°' bk \ normal(d 0 Ak — doAb, &0 2 ) ifb^K 

Since for each treatment k relative to A the impact 
of the covariate is the same, f> x Xj will cancel out 
for the comparison of treatment k relative to b when 

b^A 

normal ({d 0Ak + P x Xj) - (d 0Ab + ^Xy) , <r 0 2 ) 
= normal (dobk, Co 2 ) 

Figure 1 illustrates the results of a network meta- 
analysis of AB and AC studies according to Models 1- 
3 assuming the hazard over time follows a Weibull 
distribution, i.e. p = 0. The AB and AC planes reflect 
the pooled effect estimates of treatment B and C 
relative to A as a function of time and a covariate 
value. The vertical axis represents model parameters 
do A k with k = B } C, corresponding to the log-HR 
of treatment B and C relative to A at time t = 0 and 
covariate value X=0. The slope of the AB and AC 
planes as a function of ln(time) represents parameter 
d 1A k, the impact of time on the log HR of treatment k 
relative to A. The slope of the AB and AC planes as a 
function of the covariate value represents fi X Ah which 
is the impact of covariate X on the treatment effect 
parameter d 0A k (the scale). Figure 1A represents Model 
1 where it is assumed that the covariate is not an ef- 
fect modifier of d 0Ak and therefore (> xA k = 0- In 
Figure IB the effect of the covariate for the AB com- 
parison is different from the AC comparison (Model 
2). Figure 1C reflects Model 3 with the same effect of 
the covariate for the AB and the AC comparison. 

The random effects Models 1-3 treat multiple-arm 
trials (>2 treatments) without taking account of the cor- 
relations between the trial specific Ss that they estimate. 
Bayesian random effects models with a heterogeneity 
parameter for d 0Ak can be easily extended to fit trials 
with 3 or more treatment arms by decomposing a mul- 
tivariate normal distribution as a series of conditional 
univariate distributions [9,13]: 
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then the conditional univariate distributions for arm i 
given the previous 1,. . arms are: 
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Higher dimensional models with heterogeneity 
and covariate effects acting on multiple treatment 
effect parameters 

Models 2 and 3 have treatment effects on a scale and 
one shape parameter, and make the simplifying assump- 
tion that the heterogeneity and covariate only interact 
with the treatment effects in terms of scale. Two- 
dimensional random effects models with additional het- 
erogeneity parameters for treatment in terms of shape 
have the flexibility to capture between study variance 
regarding changes in the log HRs over time [8]. Building 
upon the higher-dimensional network meta-analysis 
models proposed by Jansen [8], a general formulation is 
a model with heterogeneity and covariate effects that act 
on 3 treatment effect parameters (i.e. 1 scale and 2 shape 
parameters) as presented at the end of this page. 

(y Ak reflect the impact of study level covariate Xj 

\fLk) 

on the pooled treatment effects in terms of scale, 8 0 jbh 
and shape Sy^ and S 2 jbh with treatment J< relative to 
control treatment b. 

Z is the between study covariance matrix with a 0 , o~ 1} o 2 
representing the heterogeneity in treatments effects £q/M> 
Sy bk and S 2 jbk respectively. p 01 , p 02 and p 12 are the 



correlations between these parameters. Although such 
a general model is very flexible to explore heterogeneity 
and inconsistency, identifiability is expected to be a 
challenge. 



Illustrative example 

An example of the models is presented for a network 
meta-analysis of treatment of advanced (stage IIIc or IV) 
melanoma. Although the interventions compared in the 
network in this example do not reflect the latest treat- 
ment options for melanoma, it provides a useful illustra- 
tion of the survival models with covariates. 

In the early stages of melanoma, surgery presents a 
potential curative option. However, for unresectable late 
stages, the mainstay of treatment is (experimental) sys- 
temic therapy. A literature search identified 10 RCTs 
evaluating the following 4 treatment groups: dacarbazine 
(DTIC) monotherapy, DTIC + Interferon (IFN), DTIC + 
Non-IFN, and Non-DTIC [15-24]. The network of RCTs 
is presented in Figure 2. 

For each treatment arm in each study the reported 
Kaplan-Meier curves were digitized (Digitizelt; http:// 
www.digitizeit.de/) and the data extracted from these 
trials were included in the network meta-analysis. In 
Figure 3 the scanned survival proportions are presented 
according to each direct comparison available. This 
aggregate data was analyzed using two-dimensional net- 
work meta-analysis models. Whilst network meta-analysis 
can be performed with a frequentist or a Bayesian 
approach, for this manuscript we focus on the Bayesian 
approach. Within the Bayesian framework, analyses con- 
sist of data, likelihood, parameters, and a model. 

The survival curves in Figure 3 can be divided into q 
consecutive intervals over the follow-up period: [t 1} t 2 ], 
(t 2 , t 3 ], . . ., (t q , t q+1 ] with ^=0. For each time interval 
m =1,2,3,. . .,.q extracted survival proportions were used 
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AC comparison 



In(time) 



Model 1: Two-dimensional treatment effect without 
covariates 



d 0AB > 0, 



ln(hazard ratio) 



d lAB >0 
d 1AC = 0 




AB comparison 



AC comparison 
In(time) 



covariate 
value 

Model 2: Two-dimensional treatment effect with 
treatment-specific covariate interactions 



ln(hazard ratio) 



d 0AB > 0 > d \AB > 0 
d 0AC > 0 > d XAC = 0 




AB comparison 



AC comparison 



In(time) 



Model 3: Two-dimensional treatment effect with constant 
treatment covariate interactions 

Figure 1 Graphical representation of two-dimensional network 
meta-analysis model for survival data without treatment-by- 
covariate interaction (A); two-dimensional network meta- 
analysis model with treatment specific covariate interactions 
(B); two-dimensional network meta-analysis model with 
constant covariate interaction (C). (See Model 1-3 and text for 
explanation). 



to calculate the patients at risk at the beginning of that 
interval and incident number of deaths. (A more specific 
explanation is provided in the Additional file 1 of this 
paper.) A binomial likelihood distribution of the incident 
events for every interval can be described according to: 

r jk rbin{p jkt ,n jkt ) 



Where r jkt 



is the observed number of events in the m th 
interval ending at time point t m+1 for treatment k in 
study /. Hjto is the number of subjects at risk just before 
the start of that interval adjusted for the subjects cen- 
sored in the interval. Pj kt is the corresponding underlying 
event probability. When the time intervals are relatively 
short, the hazard rate hj kt at time point t for treatment 
k in study / can be assumed to be constant for any time 
point within the corresponding m th time interval. The 
hazard rate corresponding to p jkt for the m th interval can 
be standardized by the unit of time used for the analysis 
(e.g. months) according to: hj kt = -ln(l -pj kt )/Atj kt where 
Atjto is the length of the interval. For the model estima- 
tion we assign this underlying hazard to time point t m+1 . 

Study date can be considered a proxy for potential 
changes over time in (background) medical care and 
other contextual factors that influence the treatment 
effects. For the included RCTs, there is variation within 
and between different types of comparisons regarding 
study date. The impact of study date on treatment 
effects was evaluated with the following meta-regression 
models: 1) random and fixed effects Weibull (p = 0) 
survival models without treatment-by-covariate inter- 
action (Model 1; Figure 1A); 2) random and fixed effects 
Weibull survival models with treatment specific covari- 
ate interaction terms (Model 2; Figure IB); 3) random 
and fixed effects Weibull survival models with a con- 
stant treatment-by-covariate interaction (Model 3; 
Figure 1C). Only Weibull models were used because the 
In (-ln(Survival)) versus ln(time) showed linear relations 
for the different studies, indicating Weibull models are 
appropriate [25]. 

Here below the prior distributions are presented as 
used for the meta-analysis. 



( normal 

\ u\Ak ) 

f$ xAk ~normal(0, 10 4 
Go-uniform^, 2) 



((o 



10 4 0 

0 10 4 

10 4 0 

0 10 4 



With constant covariate interaction, fi xAk ~ normal 
(0,10 4 ) will be replaced with /3 X ~ normal(0,10 4 ). With 
the fixed effects model, it is not necessary to define a 
prior distribution for a 0 . 
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Young et al., 2001 
Falksonetal., 1998 
Falkson etal., 1991 
Bajetta et al., 1994 
Thomson et al., 1993 




Falksonetal., 1998 



DTIC 




Falksonet al., 1998 






Cocconi et al., 1992 
Chiarion et al., 1992 
Chapman etal., 1992 




Avril et al., 2004 
Middletonet al., 2000 


Non-DTIC 





DTIC+ 
non-IFN 



Figure 2 Network of randomized controlled trials. 
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90% 



-■-Avril 2004-DTIC 
-•-Avril 2004-nonDTIC 
— -Middleton 2000-DTIC 
Middleton 2000-nonDTIC 



-Falkson 1998-DTIC+IFN 
-Falkson 1998-DTIC+nonlFN 



Time (months) 

Figure 3 Survival as observed in individual studies with different treatment comparisons. 
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The parameters of the different models were estimated 
using a Markov Chain Monte Carlo (MCMC) method as 
implemented in the WinBUGS software package [26]. 
(See Additional file 1 for the code) The first 80,000 itera- 
tions from the WinBUGS sampler were discarded as 
'burn-in' and the inferences were based on an additional 
50,000 iterations using two chains. Convergence of the 
chains was confirmed by the Gelman-Rubin statistic. 

The deviance information criterion (DIC) was used 
to compare the goodness-of-fit of the different mod- 
els [27,28]. DIC provides a measure of model fit 
that penalizes model complexity according to DIC = 
D+pD,pD = D — D [28]. D is the posterior mean 
residual deviance [28], pD is the 'effective number of 
parameters' and D is the deviance evaluated at the pos- 
terior mean of the model parameters. In general, a more 
complex model will result in a better fit to the data, 
demonstrating a smaller residual deviance. The model 
with the lowest DIC is the model providing the 'best' fit 
to the data adjusted for the number of parameters. 

Results 

Illustrative example 

In Table 1 the model fit statistics for the different mod- 
els are presented. The random effects model without 
covariates was associated with a smaller residual devi- 
ance than the fixed effects model without covariates. 
Taking into account the increased model complexity of 
the random effects approach, the DIC with the random 
effects model is also lower. We prefer the random effects 
model over the fixed effects model. (We take the pos- 
ition that in principle a random effects model is pre- 
ferred over a fixed effects model because there is often 
cause to suspect heterogeneity. Hence, in a situation 
when the DIC with the more complex random effects 
model is comparable to the fixed effects model, suggest- 
ing there is no strong evidence against the fixed effects 
model on statistical grounds, we might still prefer to use 
the random effects model to obtain a conservative 



estimate. In a situation when there is not sufficient data 
to estimate a heterogeneity parameter, a fixed effects 
model is preferred.) Adding treatment specific covariate 
interaction terms to the random effects model improved 
the residual deviance, but resulted in a similar DIC 
(1492.5) due to the greater number of model parameters. 
Treatment constant interaction effects for the random 
effects model resulted in a similar residual deviance, but 
the DIC was arguably somewhat more favorable (1490.7) 
because of the less complex model formulation. For the 
fixed effects model, extension with treatment specific 
covariate interaction terms improved both the residual 
deviance and the DIC, suggesting that factors associated 
with study date (partly) explain the heterogeneity in the 
scale related treatment effects. A fixed effect model with 
a constant covariate interaction term showed a compar- 
able model fit. 

Table 2 provides parameter estimates for the fixed and 
random effects model without covariates (Model 1), as 
well as parameter estimates for the random effects 
model with treatment specific covariate interaction 
terms (Model 2) and the random effects model with a 
constant covariate interaction term (Model 3). Although 
credible intervals for the treatment effects in terms of 
shape include the null, the point estimates (on a log HR 
scale) for DTIC + IFN and non-DTIC are sufficiently 
large to defend the two-dimensional models. (Fur- 
thermore, ignoring the treatment effects in terms of 
shape, i.e. assuming a constant HR, implies that the 
uncertainty in model shape is not captured with the 
meta-analysis models.) Adding the treatment-by-covariate 
interaction term notably affected the treatment effect 
in terms of scale (d 0Ak ). The treatment-by-covariate inter- 
action for non-DTIC vs. DTC (-0.02) was different than 
the interaction term obtained with a model with 
a constant interaction term (0.05), which implies that 
the assumption of a constant covariate interaction can 
be challenged. 

Based on the pooled treatment effects regarding the 
scale and shape of each intervention relative to DTIC 



Table 1 Goodness of fit estimates for fixed effects and random effects network meta-analysis models 



Model 


Dbar 


Dhat 


pD 


DIC 


Random effects Weibull model without covariates 


1462.4 


1432.3 


30.1 


1492.5 


Fixed effects Weibull model without covariates 


1468.9 


1443.0 


25.9 


1494.8 


Models with study data as covariate 










Random effects Weibull model with treatment specific covariate interactions 


1460.5 


1428.5 


32.0 


1492.5 


Fixed effects Weibull model with treatment specific covariate interactions 


1464.7 


1436.7 


28.0 


1492.7 


Random effects Weibull model with constant treatment covariate interactions 


1459.7 


1428.7 


31.0 


1490.7 


Fixed effects Weibull model with constant treatment covariate interactions 


1466.7 


1439.9 


26.8 


1493.5 



Dbar = posterior mean residual deviance. 

Dhat = deviance evaluated at the posterior mean of the model parameters. 
pD = the effective number of parameters as a measure of model complexity. 
DIC = Deviance Information Criteria. 



Table 2 Model parameters for fixed and random effects two-dimensional (Weibull) network meta-analysis models with and without treatment-by-covariate 
interaction 





Fixed effects model 


Random effects model (model 1) 


Random effects model with 
treatment specific covariate 
interaction* (model 2) 


Random effects model with 
constant treatment-by-covariate 
interaction* (model 3) 


Median 
of posterior 
distribution 


95% Credible 
Interval 


Median 
of posterior 
distribution 


95% Credible 
Interval 


Median 
of posterior 
distribution 


95% Credible 
Interval 


Median 
of posterior 
distribution 


95% Credible 
Interval 


Pooled estimate for difference in scale (3 0 


















DTIC + IFN vs. DTIC {d 0AB ) 


-0.16 


(-0.63; 0.33) 


-0.22 


(-0.76; 0.23) 


-0.04 


(-0.54; 0.47) 


-0.18 


(-0.58; 0.39) 


DTIC + non-IFN vs. DTIC (d 0AC ) 


-0.07 


(-0.46; 0.34) 


-0.19 


(-0.65; 0.30) 


-0.16 


(-0.66; 0.32) 


-0.10 


(-0.51; 0.35) 


non-DTIC vs. DTIC {d 0AD ) 


-0.27 


(-0.63; 0.13) 


-0.30 


(-0.88; 0.24) 


-0.17 


(-1.40; 1.11) 


-0.43 


(-1.12; 0.06) 


Pooled estimate for difference in shape £>, 


















DTIC + IFN vs. DTIC (d 1AB ) 


0.13 


(-0.08; 0.33) 


0.14 


(-0.04; 0.34) 


0.12 


(-0.07; 0.30) 


0.16 


(-0.04; 0.30) 


DTIC + non-IFN vs. DTIC (d 1AC ) 


-0.06 


(-0.23; 0.11) 


-0.02 


(-0.19; 0.15) 


-0.04 


(-0.21; 0.13) 


-0.05 


(-0.19; 0.10) 


non-DTIC vs. DTIC (d 1AD ) 


0.04 


(-0.1 7; 0.23) 


0.06 


(-0.15; 0.27) 


0.09 


(-0.11; 0.29) 


0.03 


(-0.16; 0.21) 


Estimate for covariate effect (fix) 


















DTIC + IFN vs. DTIC (fix AB ) 










0.06 


(-0.02; 0.16) 


0.05 


(-0.01; 0.12) 


DTIC + non-IFN vs. DTIC {fix AC ) 










0.06 


(-0.05; 0.18) 


0.05 


(-0.01; 0.12) 


non-DTIC vs. DTIC {fix AD ) 










-0.02 


(-0.21; 0.19) 


0.05 


(-0.01; 0.12) 


Between study variance (heterogeneity) in scale 






0.21 


(0.02; 0.56) 


0.19 


(0.01; 0.52) 


0.18 


(0.01; 0.55) 



* Covariate value set at the average across studies. 
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Figure 4 (See legend on next page.) 
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(See figure on previous page.) 

Figure 4 Treatment effects of different interventions relative to DTIC as obtained with random effects Weibull network meta-analysis 
expressed as hazard ratio (and 95% credible interval). A) Results without adjustment for differences in study date; B) Results after adjustment 
for differences in study date assuming treatment specific covariate interactions. Hazard ratios are presented corresponding to the average 
covariate value; C) Results after adjustment for differences in study date assuming constant covariate interaction. Hazard ratios are presented 
corresponding to the average covariate value. 



(d 0 Ah an d diAh with k = B,C,D, corresponding to respect- 
ively DTIC + IFN, DTIC + non-IFN and non-DTIC) the 
resultant HRs as a function of time were obtained 
according to \n(HR Ak ) = d 0 A£ + d X Ak' m M- Figure 4A 
reflects the HRs over time (along with 95% credible 
intervals) with a random effects model without adjust- 
ment for differences in study date across studies and 
comparisons. Figure 4B shows the HR over time after ad- 
justment for differences in study date using a random ef- 
fect model with treatment specific covariate interactions 
(model 2). The HRs over time are presented for the aver- 
age study date of all studies. Figure 4C shows the HRs 
over time after adjustment for study date using a random 
effect model with a constant treatment-covariate inter- 
action (model 3). Comparing Figure 4A with Figure 4B 
and 4c illustrates the effect of ignoring the variation in 
study date across the different comparisons. 

By applying the treatment effects on scale and shape 
as obtained with model 1, 2 and 3 (Table 2) to an aver- 
age scale and shape of the studies with DTIC, an 
expected scale and shape was obtained for the other 
interventions. The corresponding survival functions for 
each of the four interventions are presented in Figure 5. 
This figure allows for comparisons of survival propor- 
tions at different time points (including median survival), 
as well as comparisons of expected (i.e. mean) survival, 
which is useful for cost-effectiveness evaluations. In 
Table 3 differences in expected survival are presented for 
the different random effects models. 

Given the Bayesian approach, the probability that a 
certain treatment shows the greatest survival at different 



time points was presented based on the posterior distri- 
bution of the estimated survival proportions over time 
(Figure 6). 

To illustrate the relevance of the covariate publication 
date for this analysis, the study specific differences in 
scale of DTIC + IFN, DTIC + non-IFN, and non-DTIC 
relative to DTIC are presented as a function of publica- 
tion year in Figure 7. (These study specific estimates 
were obtained with a model similar to Model 1, but now 
assuming independent study specific differences in scale, 
rather than exchangeable effects as obtained the random 
effects model.) From Figure 7 it can be inferred that for 
the comparisons DTIC + IFN versus DTIC and DTIC + 
non-IFN versus DTIC the treatment effects in terms of 
scale show an increasing trend over the years. Further- 
more, the non-DTIC versus DTIC studies were per- 
formed at later point in time than the other studies. 
Given this imbalance, adjustment for publication date 
as a proxy for changing medical care is justified for this 
indirect comparison. 

Discussion 

In this paper, the network meta-analysis models pro- 
posed by Ouwens et al. and Jansen for survival data 
are extended with treatment-by-covariate interactions to 
explain heterogeneity and possibly to adjust for con- 
founding bias due to inconsistency [7,8]. The primary 
advantage of these evidence synthesis models for sur- 
vival data is that they do not rely on a proportional 
hazards assumption across studies and comparisons, and 
any inconsistency due to imbalance in known and 



Table 3 Difference in expected survival (in months) between interventions as obtained with random effects network 
meta-analysis model with and without treatment-by-covariate interaction 

Random effects model without Random effects model with Random effects model with 

covariate interaction (model 1) treatment specific covariate constant treatment-by-covariate 









interaction (model 2) 


interaction (model 3) 






Median 


95% 


Median 


95% 


Median 


95% 




of posterior 


Credible 


of posterior 


Credible 


of posterior 


Credible 




distribution 


Interval 


distribution 


Interval 


distribution 


Interval 


DTIC + IFN vs. DTIC 


-1.12 


(-4.20; 3.49) 


-2.46 


(-5.72; 1.91) 


-1.90 


-5.10; 2.22) 


DTIC + non-IFN vs. DTIC 


3.63 


(-1.72; 10.75) 


3.77 


(-1.04; 10.72) 


3.21 


-2.39; 9.65) 


non-DTIC vs. DTIC 


2.66 


(-2.38; 13.49) 


1.16 


(-7.51; 25.76) 


6.60 


-0.78; 21.39) 


DTIC + non-IFN vs. DTIC + IFN 


4.71 


(-1.38; 11.37) 


6.26 


(0.26; 13.29) 


5.13 


-1.15; 11.50) 


non-DTIC vs. DTIC + IFN 


3.81 


(-2.40; 14.11) 


3.46 


(-6.00; 27.88) 


8.48 


-0.65; 24.00) 


non-DTIC vs. DTIC + non-IFN 


-0.84 


(-9.15; 10.3) 


-2.79 


(-13.90; 22.17) 


3.27 


-6.72; 20.74) 
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(See figure on previous page.) 

Figure 5 Survival (and 95% credible interval) with different interventions as obtained with random effects Weibull network meta- 
analysis. A) Results without adjustment for differences in study date; B) Results after adjustment for differences in study date assuming treatment 
specific covariate interaction. Results are presented corresponding to the average covariate value across studies; C) Results after adjustment for 
differences in study date assuming constant covariate interaction. Results are presented corresponding to the average covariate value. 

V J 



measured treatment effect modifiers can be (partly) 
adjusted for. Although models to incorporate covariates 
in network meta-analysis have been presented before, 
this is the first application where treatment effects act 
on two separate parameters [10-14]. 

Since the treatment effect acts on both scale and shape 
with these network meta-analysis models, differences in 
treatment effect modifiers across studies can (in principle) 
cause heterogeneity and inconsistency in terms of both 
the scale and shape. Accordingly, treatment-by-covariate 
interactions to explain heterogeneity or minimize incon- 
sistency can be of a multidimensional nature as well. 
However, identifiability of models with treatment-by- 
covariate interactions for both scale and shape may be 
very challenging. Under the assumption that heterogen- 
eity of treatment effects only act on the scale parameter, 
we need at least 2 studies for at least one of the treatment 
comparisons in the network to estimate a corresponding 
random effects model. If also covariate interactions are 
assumed for this effect, even more studies are needed. If 
we want to estimate heterogeneity and/or treatment-by- 
covariate interactions in terms of shape, we run in the 
additional challenge of decreasing sample size and event 
counts for longer follow-up time points, thereby making 
the estimation of any shape related parameter beyond 
a fixed treatment effect challenging. To ensure model 
identifiability, it is recommended to use network meta- 
analysis models with fixed treatment effects in terms of 
shape and only heterogeneity and treatment-by-covariate 
interactions regarding treatment effects in terms of the 
scale. Given the presence of a time related treatment 
effect in these models (i.e. shape), it is unlikely that study 
characteristics, patient characteristics, and contextual fac- 
tors (which are most likely unaffected by time) have an 
impact on heterogeneity and inconsistency beyond the 
constant component of the treatment effect (i.e. the 
scale). More complex models with heterogeneity and cov- 
ariate interactions in terms scale and shape parameters 
are only expected to be feasible if a large number of stud- 
ies with sufficient number of events for longer follow-up 
times are available. 

To estimate treatment specific covariate interaction 
terms (Model 2) we need sufficient covariate variation 
across the studies for each intervention k relative to an 
overall anchor treatment (A). For models with a con- 
stant covariate interaction term (Model 3) we have fewer 
parameters to estimate and these models are therefore 
easier to identify. Furthermore, for these models we 



do not necessarily need spread in the covariate across 
studies for each type of comparison relative to A; we 
only need sufficient variation across studies comparing 
any intervention relative to treatment A. An alternative 
approach to estimate treatment specific covariate inter- 
action terms is a model with exchangeable interaction 
effects described by a normal distribution [10]. Such a 
model allows interaction estimates that shrink towards a 
common mean, thereby improving parameter estimation. 

This network meta-analysis was based on survival pro- 
portions extracted from published Kaplan-Meier curves, 
used to calculate the incident number of deaths and 
patients at risk per interval according to an algorithm 
described by Jansen (See Additional file 1 of this paper) 
[7,8]. Guyot et al. recently proposed an improved algo- 
rithm to reconstruct the data from published Kaplan- 
Meier curves, which provides a closer approximation 
of the censoring times, thereby possibly improving the 
accuracy in terms of the number of patients at risk and 
allowing greater accuracy of the uncertainty in model 
parameter estimates [29]. 

In the present paper, covariate adjustment was based on 
study-level or aggregate level data. A challenge with 
meta-regression models using study level data is that the 
association between a patient characteristic and the treat- 
ment effect of the studied interventions at the study level 
may not reflect the individual level effect-modification 
of that covariate [30,31]. Hence, the models can only be 
used with study level data when there is between- study 
or between-comparison variation with limited variation 
in effect modifiers within studies. This is typically the 
case for study design characteristics or characteristics of 
the intervention such as dose. However, in the case of 
imbalances in baseline patient characteristics across com- 
parisons, variation in these characteristics is often present 
within studies, in which case patient level data is required 
to ensure valid adjustment for inconsistency. Further- 
more, patient-level data provides a greater opportunity to 
explore differences in effects among subgroups. However, 
obtaining patient-level data for all RCTs in the network 
may be considered unrealistic. If there is only individual 
patient data available for a subset of trials, combining 
individual patient data for this subset with study level data 
from the other studies provides a useful alternative [31]. 

The proposed models can also be used to adjust for 
differences in terms of baseline risk or placebo response 
across the trials. The placebo response reflects the 
impact of all study or patient characteristics that have 
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Figure 6 (See legend on next page.) 
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(See figure on previous page.) 

Figure 6 Probability of greatest survival up to a certain time point (AUC) with different interventions as obtained with random 
effects Weibull network meta-analysis. A) Results without adjustment for differences in study date; B) Results after adjustment for differences 
in study date assuming treatment specific covariate interaction. Results are presented corresponding to the average covariate value across 
studies; C) Results after adjustment for differences in study date assuming constant covariate interaction. Results are presented corresponding 
to the average covariate value. 



an impact on the outcome; the placebo response sum- 
marizes the impact of prognostic factors (i.e. the study 
effect). If there is an association between the placebo 
response and treatment effects, it implies that some or 
all of the study and patient characteristics that are prog- 
nostic factors of the outcomes are also treatment effect 
modifiers. Adjustment for the placebo response partly 
mitigates inconsistency or bias due to an imbalance in 
effect modifiers across comparisons [9,10]. This may be 
important given the challenge of adjusting for multiple 
differences in effect modifiers using study-level data and 
the limited feasibility of accessing individual patient data 
for all trials. A limitation of adjusting for baseline risk in 
a network meta-analysis is the (theoretical) possibility of 
introducing collider stratification bias [9]. 

Although network meta-analysis where the treatment 
effects contain a time related aspect has obvious advan- 
tages in terms of model fit to the data, the presenta- 
tion of the results might need some familiarization. 
The advantage of presenting HRs as a function of time 
(Figure 4) is that time varying treatment effects can be 
identified, but a possible disadvantage is the challenge to 
identify whether one treatment is overall favourable over 
another, especially when the HR curves cross. Pooled 
treatment specific survival curves as shown in Figure 5 



provide information on commonly understood concepts 
like median survival and survival at different time points. 
A disadvantage, however is that a baseline survival curve 
with treatment A needs to be defined in order to trans- 
late the HR curves as obtained with the network meta- 
analysis into (pooled) survival curves by treatment. 
In the current example the pooled curve with treatment 
A was based on the average of study specific nuisance 
parameters for scale and shape of all treatment A (i.e. 
DTIC) controlled trials. 

Within a Bayesian framework it is possible to calculate 
the probability of being the best treatment out of all 
those treatments assessed, or second best, third best, etc. 
[32]. With the time varying treatment effects obtained 
from the network meta-analysis models presented in this 
paper there are different options to create probabilistic 
summaries of treatment effects. We can either present 
a probabilistic summary based on a collapsed measure 
of the time varying treatment effects, or present prob- 
abilistic summaries as a function of time. Examples of 
the first category are based on the median survival or 
expected survival obtained from the pooled survival 
curves. Probabilistic summaries over time can be based 
on the HR at each time point, the survival proportion at 
each time point, or the area under the survival curve 
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Figure 7 Independent study specific differences in scale of DTIC + IFN, DTIC + non-IFN, and non-DTIC relative to DTIC as a function of 
publication year. The associations between treatment effects in terms of scale and publication year as obtained with model 2 are presented 
with the dashed lines. 
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(AUC) to the left of each time point. Figure 6 is based on 
this concept and summarizes the cumulative treatment 
effect. Additional research is required to understand 
the sensitivity of these different probabilistic measures to 
the estimated underlying basic model parameters for the 
treatment effect, and provide some guidance. 

A multidimensional network meta-analysis model for 
survival data is extremely useful for cost-effectiveness 
analysis because estimates of expected survival differ- 
ences of competing interventions are less likely to be 
biased [7,8]. The extension of the two-dimensional 
model to include covariates allows for an evaluation of 
patient subgroups for which the clinical or cost effective- 
ness of the technology might be expected to differ from 
the overall population. Although the current analysis 
focused on overall survival, this method can be applied 
to any time to event outcomes. It may be of interest 
to extend the model to evaluate both progression-free 
and overall survival simultaneously using a multivariate 
approach, as has been proposed by Welton et al. [33]. 
This will avoid any inconsistency between clinical evi- 
dence synthesis and economic evaluations based on 
models with differences in quality of life before and after 
disease progression. 

Conclusion 

Adding treatment-by-covariate interactions to multidi- 
mensional network meta-analysis models for published 
survival curves can be worthwhile to explain systematic 
differences across studies and to reduce inconsistencies. 
An additional advantage is that heterogeneity in survival 
data can be addressed. These models are not only useful 
for comparative effectiveness evaluation, but also pro- 
vide an opportunity to ensure consistency within a cost- 
effectiveness analysis. In the Additional file 1 the data 
requirements to perform analysis with these kinds of 
models are outlined. 

Additional file 



Additional file 1: Data requirements to perform network 
meta-analysis of published survival curves. 
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